#include "stdafx.h"

/*
 * -- SuperLU routine (version 2.0) --
 * Lawrence Berkeley National Lab, Univ. of California Berkeley,
 * and Xerox Palo Alto Research Center.
 * September 10, 2007
 *
 */
/*
 * File name:	dlangs.c
 * History:     Modified from lapack routine DLANGE
 */
#include <math.h>
#include "hnum_pdsp_defs.h"

namespace harlinn
{
    namespace numerics
    {
        namespace SuperLU
        {
                double dlangs(char *norm, SuperMatrix *A)
                {
                /* 
                    Purpose   
                    =======   

                    DLANGS returns the value of the one norm, or the Frobenius norm, or 
                    the infinity norm, or the element of largest absolute value of a 
                    real matrix A.   

                    Description   
                    ===========   

                    DLANGS returns the value   

                       DLANGS = ( max(abs(A(i,j))), NORM = 'M' or 'm'   
                                (   
                                ( norm1(A),         NORM = '1', 'O' or 'o'   
                                (   
                                ( normI(A),         NORM = 'I' or 'i'   
                                (   
                                ( normF(A),         NORM = 'F', 'f', 'E' or 'e'   

                    where  norm1  denotes the  one norm of a matrix (maximum column sum), 
                    normI  denotes the  infinity norm  of a matrix  (maximum row sum) and 
                    normF  denotes the  Frobenius norm of a matrix (square root of sum of 
                    squares).  Note that  max(abs(A(i,j)))  is not a  matrix norm.   

                    Arguments   
                    =========   

                    NORM    (input) CHARACTER*1   
                            Specifies the value to be returned in DLANGE as described above.   
                    A       (input) SuperMatrix*
                            The M by N sparse matrix A. 

                   ===================================================================== 
                */
    
                    /* Local variables */
                    NCformat *Astore;
                    double   *Aval;
                    int      i, j, irow;
                    double   value, sum;
                    double   *rwork;

                    Astore = (NCformat*)A->Store;
                    Aval   = (double*)Astore->nzval;
    
                    if ( SUPERLU_MIN(A->nrow, A->ncol) == 0) {
	                value = 0.;
	
                    } else if (lsame_(norm, "M")) {
	                /* Find max(abs(A(i,j))). */
	                value = 0.;
	                for (j = 0; j < A->ncol; ++j)
	                    for (i = Astore->colptr[j]; i < Astore->colptr[j+1]; i++)
		                value = SUPERLU_MAX( value, fabs( Aval[i]) );
	
                    } else if (lsame_(norm, "O") || *(unsigned char *)norm == '1') {
	                /* Find norm1(A). */
	                value = 0.;
	                for (j = 0; j < A->ncol; ++j) {
	                    sum = 0.;
	                    for (i = Astore->colptr[j]; i < Astore->colptr[j+1]; i++) 
		                sum += fabs(Aval[i]);
	                    value = SUPERLU_MAX(value,sum);
	                }
	
                    } else if (lsame_(norm, "I")) {
	                /* Find normI(A). */
	                if ( !(rwork = (double *) SUPERLU_MALLOC((size_t) A->nrow * sizeof(double))) )
	                    SUPERLU_ABORT("SUPERLU_MALLOC fails for rwork.");
	                for (i = 0; i < A->nrow; ++i) rwork[i] = 0.;
	                for (j = 0; j < A->ncol; ++j)
	                    for (i = Astore->colptr[j]; i < Astore->colptr[j+1]; i++) {
		                irow = Astore->rowind[i];
		                rwork[irow] += fabs(Aval[i]);
	                    }
	                value = 0.;
	                for (i = 0; i < A->nrow; ++i)
	                    value = SUPERLU_MAX(value, rwork[i]);
	
	                SUPERLU_FREE (rwork);
	
                    } else if (lsame_(norm, "F") || lsame_(norm, "E")) {
	                /* Find normF(A). */
	                SUPERLU_ABORT("Not implemented.");
                    } else
	                SUPERLU_ABORT("Illegal norm specified.");

                    return (value);

                } /* dlangs */

        };
    };
};